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We numerically demonstrate the formation of carrier field shocks in various dispersive media for 
a wide variety of input conditions using two different electric field propagation models. In addition, 
an investigation of the impact of numerous physical effects on carrier wave shock is performed. It is 
shown that in many cases a field shock is essentially unavoidable and therefore extremely important 
in the propagation of intense long wavelength pulses in weakly dispersive nonlinear media such 
as noble gases, air, and single-crystal diamond. The results presented here are expected to have 
a significant impact in the field of ultrashort nonlinear optics, attosecond pulse generation, and 
wavepacket synthesis where the use of mid-IR wavelengths is becoming increasingly more important. 


I. INTRODUCTION 

Traditionally, the field of ultrafast nonlinear optics uses 
near IR lasers for the study of the interaction between 
atoms and molecules with the intense electric field of the 
pulse, i.e. 800 nm for Titanium Sapphire systems utiliz¬ 
ing the chirped pulse amplification technique n. Typi- 
cally the duration of such wavepackets is on the femtosec¬ 
ond time scale containing multiple electric field cycles. 
Therefore, the nonlinear propagation dynamics can be 
successfully modeled using Schrodinger type equations, 
the most well known of which is the nonlinear envelope 
equation (NEE) [2j [22]. While NEE based models are 
well suited to describe most effects relevant to the electric 
field envelope, they fail by definition to describe carrier 
wave related effects, most notably for the present work, 
carrier wave steepening. 

Carrier wave self-steepening and field shock formation 
have been shown to exist in an ideal dispersionless media 
through numerical integration of Maxwell’s equations [3|, 
and have very recently been predicted in realistic disper¬ 
sive media at mid and far-IR wavelengths [4] using the 
modified Kadomtsev-Petviashvili equation (MKP) [IHS]. 
The key factor for this latest result was the use of longer 
wavelengths that are beneficial both in terms of a higher 
ionization threshold and an almost flat dispersion land¬ 
scape. The resulting propagation regime is strongly non¬ 
linear and weakly dispersive, which proves to be ideal for 
the observation of extreme carrier shocks. On the other 
hand, the high dispersion and weak nonlinearity found in 
most transparent media at shorter wavelengths prevents 
carrier wave steepening from manifesting. 

As was shown in Ref. [4], the canonical description 
for pulse propagation in the strong nonlinear - weakly 
dispersive regime is the carrier wave resolved modified 
Kadomtsev-Petviashvili (MKP) equation. MKP models 
propagate an electric field rather than a wavepacket en¬ 
velope, hence they are related to the unidirectional pulse 
propagation equations (UPPE) [9] which are widely used 


in computational nonlinear optics as an alternative to en¬ 
velope based models HO]. In addition, the MKP model 
is able to predict a variety of propagation related effects 
that reshape the electric field in spectacular ways, pro¬ 
ducing exotic ’’shark-fin” and ’’top-hat” waveforms. 

The reaiization that mid and far-IR puise propaga¬ 
tion is very different than near-IR fiiamentation has far 
reaching impiications since research in the uitrashort op¬ 
tics community has trended towards ionger waveiengths 
in recent years. In Ref. m, 2-/im puises were used to 
produce a seif-compressed carrier-enveiope phase stabie 
few-cycie puise through fiiamentation. Eurthermore, a 
recent work by Berge et. al. [12] iliustrated that seif- 
compression efficiency is improved at ionger waveiengths; 
in Berge’s study, carrier waveiengths of up to S-jam were 
simuiated in argon. Extreme spectrai broadening that is 
observed in the mid-IR regime has been shown to be abie 
to produce keV x-rays [13] and opens the way for attosec¬ 
ond m and even zeptosecond puise generation [15] , and 
wave-form synthesis m- 

Given the findings for mid and far-IR waveiength 
puises in [H [12], it is evident that a reassessment of 
a wide variety of modeis in muitipie discipiines in the 
fieid of extreme nonlinear optics is necessary. Eull elec¬ 
tric field propagators are mandatory for long wavelength, 
high nonlinearity - low dispersion pulse propagation since 
carrier field shocks are likely to occur in various media 
[4]. It is clear that propagation effects have to be taken 
into account when conducting experiments in this regime 
since the field profile can no longer be considered station¬ 
ary. Currently, however, the theoretical understanding of 
the strong nonlinear - weakly dispersive regime is limited. 
Euture studies and analysis of this area will be important 
as more mid and far-IR mJ femtosecond lasers become 
commercially available. 

In this work we present an extensive numerical study 
of the strongly nonlinear - weakly dispersive propaga¬ 
tion regime found for long wavelength (mid-IR and far- 
IR) femtosecond pulses for a variety of materials with 
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two different computational models, the unidirectional 
pulse propagation equation and the modified Kadomtsev- 
Petviashvili equation. In our numerical experiments we 
use a variety of input pulses and materials, while we also 
isolate various physical effects in order to study their im¬ 
pact on carrier wave shock formation and sustenance. 
Finally, a comparison of the two theoretical models is 
used to verify our results. 


II. NUMERICAL MODELS 


In computational nonlinear optics, solving the full 
spatio-temporal Maxwell model is known to be extremely 
expensive in computation time and hence is impractical. 
Therefore, numerical studies of ultrashort pulse propa¬ 
gation are traditionally conducted either by using enve¬ 
lope type models like the nonlinear Schrodinger equa¬ 
tion [2] , or full electric field propagators such as the uni¬ 
directional pulse propagation equation (UPPE) [9]. In 
section |II A| we discuss the most general unidirectional 
field propagation method, the UPPE system. In sec¬ 


tion IIB , we derive the generalized modified Kadomtsev- 
Petviashvili (GMKP) equation from the UPPE model, 
show Keldsyh ionization rates for various media at 8-/im, 
and analyze the relevant length scales of the MKP equa¬ 
tion. Throughout this letter, results from the GMKP 
and UPPE models will be presented and compared. 


A. Fully spectral solvers: UPPE 

In fully spectral solvers, the native representation of 
the optical field is stored in arrays representing spec¬ 
tral amplitudes that are functions of angular frequency 
cj, and transverse wavenumbers k^^ky. For example 
z) would stand for the spectral amplitude in 
the plane-wave expansion of the field at a given propa¬ 
gation distance 2 ;. 

The spectral nature of the solver makes it possible to 
solve the linear propagation problem exactly. For a given 
medium permittivity e(ci;), the linear propagation prop¬ 
erties are encoded in the propagation constants of the 
modal fields. In the case of a homogeneous medium these 
are the well-known plane waves, and the propagation con¬ 
stant depends on both the angular frequency and on the 
transverse wavenumbers: 

kx, ky) = y^w 2 e(w)/c 2 -kl-kl . ( 1 ) 

The permittivity is in general complex-valued, so that 
not only the chromatic dispersive properties, but also 
absorption can be accounted for exactly. 

Medium properties other than linear enter the prop¬ 
agation model through the constitutive relations in 
Maxwell equations. Thus, polarization P(x, 2 :, t) and 
the current density J{x^y^z^t) must be given as func¬ 
tionals of the electric field E{x^y^z^t). For the sake of 


brevity, we will not show the dependence of either P or 
J on the electric field, because from the point of view 
of the propagation solver, the nature of this dependence 
is irrelevant. Let us just note one important aspect and 
that is that the medium response is calculated from the 
history of the electric field E{x^y^z^t) independently at 
a each spatial point. 

The following is an exact system of equations that de¬ 
scribes evolution of modal amplitudes along the ^-axis for 
the forward and backward propagating field. It is a cou¬ 
pled pair of unidirectional pulse propagation equations 
(UPPE) [9]: 


s=l,2 

( 2 ) 


which describes the forward propagating wave, and 


lUJ 


UJ 




( 3 ) 


which represents the backward propagating radiation. 

The above UPPEs are mutually coupled through the 
medium response represented by the polarization and 
current density terms. As noted above, they require the 
knowledge of the full (i.e. forward and backward) elec¬ 
tric field. This is where the unidirectional approximation 
must be adopted. It requires that the medium response 
can be accurately calculated from only the forward por¬ 
tion of the field, i.e.: 

P{E_ + E+), J{E_ + E+) ^ P{E+), J{E+) (4) 


If this is the case, then it is sufficient to solve only 
the UPPE representing the forward wave. Obviously, 
in practice all pulse propagation simulations, indepen¬ 
dently of the propagation model they utilize, require this 
approximation. The physical meaning of “unidirection¬ 
ality” was discussed by Kinsler in pTl[2Q] and a practi¬ 
cal way to quantify and test its limitations is described 
in [ 21 ]. The canonical form of the scalar UPPE is as fol¬ 
lows (for brevity, we consider the current term grouped 
with the nonlinear polarization and omit the component 
and direction subscripts): 

^z^kxiky (^: ^) ~ 'Id^ikx-) ky-) Ld)Ekx,ky (^: ^) 

+ iQ{kj,, ky,uj)Pkx,ky (^, (5) 

with the exact linear propagator 

K{ky;,ky,u) = ^e{uj)/-kl-kl 


(6) 
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and the nonlinear response coupling 

^2 

Uj) = . (7) 

2eoC^^uj‘^e{uj)/c^ - - k‘^ 

This is in fact the form most often utilized in practical 
simulations. 


B. MKP 


While the canonical UPPE predicts carrier field shock 
for long-wavelength high power lasers, this effect and 
associated scales are hidden in the general form of the 
model itself. To reveal the physics in a more explicit way, 
we make several approximations to equation (§ that are 
valid for a large class of filamentation problems yet re¬ 
duce the model to a simpler and easier to analyze form. 

First, we apply the following Taylor expansions to 
K{kx,ky,uj) diiid Q{kx,ky,uj) 


kzij^i kx') ky^ k{uj^ 


k 1 kx 1 ky^ k{uj^ 


ki 


2k‘^{uj) Sk^{uj) 
+ ' 


kl 


2k‘^{(jj) 


with k^ = k‘^ y- ky^ and k{uj) = Then our 

expanded UPPE is 


dzEk^,ky (^, z) = ik{uj)Ek^^ky z) 


ick‘^\ 


2n{uj)L 


~Eka,,ky{^, + 


lUJ 


2€ocn{uj) 


Pk^,ky (^,^) 


ick‘^\ k^\ . . 

Ek,,ky {UJ,Z) 


8n(cj)cj k^(cj) 
iuj k\ 
4eocn(ci;) kP‘(uj) 


+ h.o.t. (8) 


where h.o.t. denotes the higher order terms in the re¬ 
spective expansions of kz and k~^. Restricting ourselves 
to paraxial beam propagation implies that only the first 
three right-hand-side terms must be retained. Further¬ 
more, n~^{uj) in the remaining terms can be Taylor ex¬ 
panded around a reference frequency ujr such that 


^zEkx,kyiy^s E) ~ 
ick\ 

2n{ujR)uj 

iuj 

2eocn{ujR) 


ik{uj)Ek^^ky 

_ n'{ujR){uj - Ur) 
n{ujR) 

_ n'{uR){u - Ur) 
n{uR) 


Ek^,ky 

Ek^,ky 

(9) 


Now we invoke a restriction on the medium and spectral 
bandwidth of pulses admitted to our model 


|n'(wfl)(w - Wfl)| <C n(wfl). 


This is akin to the slowly-evolving-wave approximation 
(SEWA) 

\n'{ujR)ojR\ <C n{ojR), 

which is used in envelope models to accurately describe 
light propagation down to the single-cycle regime [22]. 
Our approximation is somewhat more demanding as it 
puts additional restrictions on the medium when propa¬ 
gating high harmonic generating pulses; as the harmonic 
order increases, u will deviate more from ur and hence 
n'{uR) will need to be small enough to account for this. 
After approximations, the resulting propagation equation 
is 


dzEk,^ky{^,z) =ik{u)Ek,^ky{^,z) 


ick^\ rn / X 

-Ek^^ky i^,z) + 


^u 


2n{uR)u 


2eocn{uR) 


Pk^,ky{OJ,Z). 


( 10 ) 


The following dispersion relation approximation is uti¬ 
lized 


k{u) = au^ -h qu. 

u 


(11) 


While no explicit ur is present in the dispersion relation, 
equation 0 is a local approximation as the parameters 
a, 6, and q are chosen to accurately reflect the full disper¬ 
sion profile (e.g. Sellmeier, etc.) near a specific frequency 
Ur. Each example in this work considers a single laser 
pulse launched with a carrier frequency uq corresponding 
to a wavelength Aq. Therefore, we take the reference fre¬ 
quency as Ur = uq and let no = u^ur). Upon substitu¬ 
tion into equation (10) and transforming to time-domain 


real-space, the propagation model simplifies to 


dr (d,E - adlE + P — drP) 

\ 2eocno J 

+ bE=pA^E, (12) 

2no 

where r = t — qz and the fields are functions of (x, r, z). 
At this point it is easy to separate the nonlinear current 
density from the nonlinear polarization by replacing drP 
with drP P J- Furthermore, we split J into a lossless 
plasma defocusing term and a nonlinear absorption term 


J — Jplas T Jabsi 


^rJplas pE^ 

rrie 

Jabs _ W{I)Ui 


{Pnt P)E<, 


eocriQ I 
drP=W{I){pnt-p). 


(13) 

(14) 

(15) 

(16) 


where / = eocno < > is the intensity of electric field, 

Ui is the potential, pnt is the neutral plasma density, and 
W is the ionization rate of free charge p. We also use an 
instantaneous Kerr nonlinear response such that 


P = = ^eo(eocnon2)-B^. 
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FIG. 1. (Color online) Keldysh photoionization rates for air 
(approximated by that of oxygen), diamond, and xenon at a 
carrier wavelength of S/xm. 


Its instructive to normalize and non-dimensionalize this 
equation with 



^ r 

T = ujqt, r = —, 
Wo 


z = 


2: 

V 


where Iq is the initial intensity, wq is the beam diameter, 
and L is an arbitrary length scale. After substitutions 
we obtain 


dr d~M - 


^‘^Lcsd 


d^U + 


AL 


= 


( 22 ) 


with length scales 


This produces a propagation model with the electric field 
in units of K/m, 


Ldf = 


k{uJo)wQ 


Lnl = 


idon2lo ’ 


^CSD 


1 

12aujQ 

(23) 


v(3) 1 

dr [ d,E - ad^E + ^drE^ + 


2cno 


2eocno 


Jabs 


- p E = —A^E. (17) 

zeoTiocme J 2no 


Rescaling E by 


E = 


epcnp y/^ ^ 

2 J 


and removing tilde’s after substitution, we obtain the 
generalized modified Kadomtsev-Petviashvili (GMKP) 
optical filamentation equation 


dr (^d,E - adlE + ^drE^ + KUp, 

+ (6 + bpiasip)) E = -— A_\_E, 
Znp 


(18) 


The diffractive length scale Lqf and nonlinear length 
scale Lnl are standard length scales in nonlinear optics 
that describe the propagation of an electric field enve¬ 
lope. On the other hand, the carrier shock dispersion 
length scale Lqsd describes the propagation of the elec¬ 
tric field, not the envelope; in particular, Lcsd measures 
the distance over which carrier wave steepening dissi¬ 
pates. With the MKP dispersion relation 0 . Lcsd 
can alternatively be written as 


k{3ujp)-3k{ujp)^ 

which clearly relates the carrier shock dispersion to the 
coherence between the fundamental and the third har¬ 
monic. This important observation connects carrier wave 
shock to harmonic generation and harmonic walk-off. A 
complete discussion of Lcsd is given in section [V| 


with 


C. UPPE vs GMKP 


bplas{p) 


e 


2 


Zepupcrrie 




Afabs (P? 


W{I)Ui{pnt-p) 

21 


(19) 

( 20 ) 


At long-wavelengths, a full Keldysh ionization response 
W{I) is preferable to the multiphoton approximation 
[ 12 ]. In Fig. we show Keldysh photoionization rates 
for air (approximated by oxygen), diamond, and xenon 
at a carrier wavelength of Spm. These ionization rates 
are computed from formulas derived in [23] . 

The modified Kadomtsev-Petviashvili (MKP) equation 
is generated by dropping the nonlinear current density 
terms and ignoring anomalous dispersion (6 = 0). With 
these simplifications, we have 


dr ( d^E - adlE + '^drC ) = ( 21 ) 


3c 


2np 


The essential difference between the UPPE and GMKP 
models is in the way dispersion is modeled; the UPPE 
accurately represents dispersion over a global frequency 
range whereas the MKP dispersion, as shown in Eq. 0, 
is a local approximation whose accuracy depends on the 
the carrier frequency. In addition, the UPPE model takes 
into account vibrational Raman, the full complex x re¬ 
lation including transmittance, and plasma losses due to 
a finite electron collision time. As will be shown later 
on, even though UPPE is a more complex and com¬ 
plete model for nonlinear wave propagation, the physics 
that govern carrier shock formation can be adequately 
described by the simpler GMKP model. 

Unless stated otherwise, results presented in this work 
are obtained using the GMKP model. However, the 
UPPE is used to compare GMKP results with the com¬ 
plete model (Eig.[T7|) and to investigate effects either not 
accounted for in the GMKP model, or those that would 
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FIG. 2. Refractive index of single-crystal diamond as a func¬ 
tion of wavelength. 


FIG. 3. (Golor online) Imaginary value of y of single-crystal 
diamond as a function of wavelength. Note the two-phonon 
absorption peaks between 4 pm and 6 pm. 


push the GMKP beyond its validity limits. In particu¬ 
lar, the UPPE is used to study the impact on carrier field 
shock due to Raman (Fig. [^, high numerical aperture 
focusing (Fig.[^, and artificially high plasma generation 

(Fig.[^. 


III. MEDIA AND PULSE DESCRIPTION 

The numerical experiments are conducted in a variety 
of modeled materials ranging from atomic gases to solid 
crystal diamond. The common denominator in all cases 
is the fiat dispersion curve at mid and long IR wave¬ 
lengths (from ^ 4-pm up to 8-pm) that is essential for 
the formation of a carrier field shock. 

In our previous work we presented the first prediction 
of carrier shock formation in noble gases [4]. In this 
follow-up work, we expand the study of field shock for¬ 
mation to other materials such as air and single-crystal 
diamond. However, since noble gases are a very attrac¬ 
tive option for field shock formation, a complementary 
study was performed in order to accurately predict the 
regime in which field shocks are likely to occur (section 
V). The dispersive and nonlinear properties of the noble 
gases were taken from [ HHS ]. 

Air is a more complex gas since it consists of oxy¬ 
gen and nitrogen, and exhibits vibrational and rotational 
nonlinear responses. Nevertheless, it is another good can¬ 
didate for carrier shock formation due to the fiat dis¬ 
persive curve at longer wavelengths [26]. The nonlinear 
refractive index is n 2 = 3.2 x 10“^^ cm^/W [27], and 
the ionization potential is taken to be that of oxygen 
Ui = 9.0eV [2], which is very close to Xenon. Although 
air has a Raman response that is well known [28] [29] , it 
will not be taken into account in this work. The effect of 
Raman will be investigated for the case of single-crystal 
diamond, as will be discussed below. 

In addition to gaseous media, solids can potentially be 
suitable for producing carrier shocks due to their higher 
nonlinear responses. Concerning nonlinear field reshap¬ 


ing, diamond is in fact a close crystalline analogy to an 
atomic gas like Xenon. Single crystal diamond is a mate¬ 
rial not often used in the field of nonlinear optics, however 
it is ideal for the purpose of studying carrier field shocks. 
The main reasons for this are: First, diamond has a rela¬ 
tively high nonlinear refractive index of n 2 = 1.3 x 10“^^ 
cm^/W [30] (about one order of magnitude higher than 
fused silica glass), which is important since field shock 
formation is essentially a nonlinear phenomenon. Second, 
diamond has an almost fiat dispersion curve at longer 
wavelengths [SUES as can be seen in Fig.|^ In addition, 
the imaginary part of the susceptibility x is almost fiat if 
not for the weak two phonon absorption lines around 5- 
yum [33], shown on Fig.[^ Third, diamond has a relatively 
high ionization potential for a solid crystalline material of 
Ui = 5.5eU, which makes the ionization relatively weak 
at longer wavelengths. Diamond will also be used to as¬ 
sess the effect of vibrational Raman on field shock forma¬ 
tion. The Raman frequency uj = 2.5099 x 10^^ s“^ and 
time scale r = 4.2 ps are taken from [34l|35]. Our simu¬ 
lations have included most of the relevant effects needed 
to model nonlinear propagation of femtosecond pulses in 
diamond. These simulations support our argument that 
diamond is an attractive material for electric field steep¬ 
ening (see Figj§ for example). Furthermore, we expect 
that results presented here will motivate future exper¬ 
iments utilizing diamond as a material in the study of 
nonlinear phenomena including but not limited to field 
shocks. 

Other, less exotic solids, like fused silica or BK7 glass, 
have proven to be unsuitable to produce a field shock. 
The reason for this is the high and even anomalous dis¬ 
persion they exhibit in the IR, which will be highlighted 
later in section V. 

Throughout this work, a Gaussian spatio-temporal 
laser pulse is used, with a wavelength of Aq = 8 jam 
and a pulse duration of Tp = 40 fs at 1/e^ radius. Unless 
stated otherwise, the starting peak intensity was set to 
be Jo = 4 X 10^^ W/cm^, while the beam is collimated 
with a waist of icq = 2 mm at 1/e^ radius for air, and 
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wq = 50 jivu for diamond. All simulations, with the ex¬ 
ception of those corresponding to Fig. are done in 
three dimensions (3D): with assumption the of a cylin- 
drically symmetric solution in r, time t, and propagation 
distance z. 


IV. FIELD SHOCK FORMATION AND 
EVOLUTION 


In Fig. I^we can see the generated carrier wave shock 
in air (a, b) and diamond (c, d). The input power (beam 
waist) and propagation distances are different in each 
case to compensate for the material properties but the 
end results are very similar. In both cases (a, c) the elec¬ 
tric field of the laser pulse is undergoing extreme steep¬ 
ening and is reshaped into a ’’shark-fin” like temporal 
structure. After additional propagation, both fields con¬ 
tinue to evolve and are transformed into almost top-hats 
(b, d). These extreme nonlinear dynamics acting on the 
carrier wave, in two totally different materials, indicate 
that the mechanisms leading to carrier wave shock for¬ 
mation are universal in the strongly nonlinear - weakly 
dispersive long wavelength regime. 

It is important to stress that in this long wavelength 
- highly nonlinear - low dispersive regime, propagation 
effects must be accounted for. As we can see in Fig. 
where the evolution of the electric field of our wavepacket 
in air and diamond is shown, the field shape is in fact con¬ 
stantly evolving as the pulse is propagating in z. In both 
cases, after the initial steepening which can be seen in (a) 
at 2 : = 30 cm and (b) at 2 : = 80 /im, the field is constantly 
reshaping before finally taking on a top-hat shape at (a) 
z = 12 cm, and (b) 2 : = 200 /im. Obviously the rate 
at which the field evolves is directly proportional to the 
characteristic lengths Lcsd and The continuous 

broadening of the spectrum and walk-off of the numerous 
harmonics results in a dynamically evolving wave-form, 
which can change drastically over just a few microns of 
propagation. Therefore, any experiments and numerical 
models related to nonlinear field shaping [T3HT6] should 
be reassessed given the results presented here. 

Another interesting observation in this low ionization 
regime is that collapse can still be arrested, this time 
due to the walk-off of the self-generated higher harmon¬ 
ics. This is clearly observed in Fig. |^a) where the peak 
intensity along propagation in single-crystal diamond is 
shown (red continuous line). The extremely sharp field 
gradients which evolve due to the low dispersion - high 
nonlinearity lead to even more spectral broadening. This 
eventually results in the spreading of a significant portion 
of the total energy over a large spectral region, effectively 
dispersing the wavepacket and arresting collapse. Note 
that the role of dispersion is key; on the one hand, very 
weak dispersion allows for a carrier wave shock to man¬ 
ifest itself, on the other hand, non-zero dispersion is ar¬ 
resting the collapse and eventually breaking up the field 
profile itself. The role of plasma in this collapse arrest is 




t(fs) 



t(fs) 



FIG. 4. (Color online) Electric field of an 8-/im pulse propa¬ 
gating in air (a and b), and diamond (c and d). Black dashed 
line: initial field shape, red lines: field shape at (a) 2 : = 26.3 
cm and (b) 2 ; = 61.5 cm, for air, and (c) 2 ; = 60 fim. and (d) 
2 ; = 170 /xm, for diamond. 


actually very small. As observed in Fig. [^b) (red con¬ 
tinuous line), the generated plasma density is essentially 
kept under ^ 10^^ cm“^, which is rather weak in the 
filament at ion regime. To further test this, we switched 
plasma generation off entirely (J = 0) and observed that 
the peak intensity of the wavepacket was essentially unaf¬ 
fected (Fig.j^a)). That is, carrier shock dispersion is the 
primary mechanism by which critical collapse is arrested, 
in contrast to a shorter wavelength regime in which the 
main arrest mechanism is due to plasma. 

The role of field dispersion and envelope dispersion are 
significantly different here. When wavepacket propaga¬ 
tion is simulated using an envelope type nonlinear Kerr 





































7 



FIG. 5. On-axis electric field evolution in (a) air, and (b) 
single-crystal diamond at Ao = 8 /xm. In both cases, Iq — 
4.0 X 10^^ W/cm^ and tp — 40 fs, with wq = 2 mm for air 
and wq — 50 /xm for diamond. The electric fields initially 
evolves into a shark-fin shock with harmonic walk-off further 
reshaping the fields over longer distances. 


response (P ^ /P), instead of a carrier-resolving Kerr re¬ 
sponse (P ^ P^), the collapse arrest is caused by optical 
field ionization and plasma defocusing. This can be seen 
in Fig. where the peak intensity (a) and peak electron 
density (b) are plotted in the blue dash-dotted lines. The 
results shown in Fig. [^showcase again the importance of 
the use of field propagators for mid-IR wavelengths; en¬ 
velope type models are unable to properly describe the 
evolution of the wavepacket, predicting collapse arrest 
due optical field ionization rather than carrier shock dis¬ 
persion, and overestimating plasma generation and non¬ 
linear losses in the process. 

As shown in our previous work [4] , focusing conditions 
can play an important role in the carrier wave shock for¬ 
mation and evolution. In Fig. [^we can see the effect of 
an initial wavefront curvature on shock formation of our 
wavepacket propagating in single-crystal diamond, both 
in the temporal domain (a, c) and the spatial domain (b, 
d). Focusing the beam with a curvature of / = 50 /xm (c, 
d) will lead to earlier carrier shock formation and a top- 
hat field at z = 71 /xm. In contrast, the collimated beam 



z (mm) 



z (mm) 

FIG. 6. (Golor online) Gritical collapse of a Ao = 8 /xm Gaus¬ 
sian wavepacket propagating in single-crystal diamond with 
an input power of Pin = 1.46Pcr- Gollapse arrest due to 
plasma (P ~ IE) and carrier shock dispersion (P ~ P^) are 
illustrated in plots of (a) peak intensity vs. propagation dis¬ 
tance, and (b) peak plasma vs. propagation distance. 
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FIG. 7. (Golor online) Electric field (first column) and beam 
waist (second column) for collimated (first line) and focused 
(second line) beam in single-crystal diamond. Black line: ini¬ 
tial condition, red line: after propagation. 


















































































FIG. 8. Electric field at several radial positions after propa¬ 
gating a Ao = 8 /im collimated Gaussian beam for 170 /xm in 
single-crystal diamond. 


is still evolving at z = 200 /xm, and has not reshaped 
into a top-hat yet. The obvious explanation of this be¬ 
havior is the faster intensity increase due to strong linear 
focusing, which accelerates the harmonic generation and 
field steepening. Similar control over the field shape was 
predicted for Xenon in [4] , and is also observed in air and 
the rest of the noble gases (data not shown). Note that 
Fig. □ was generated using the UPPE model with ion¬ 
ization switched off. This was done to eliminate changes 
in the beam profile due to nonlinear losses and plasma 
defocusing. 

The transition from a collimated beam to a focused one 
has further implications in the spatial domain as well. 
As we can see in Fig. (d), the initial curvature leads 
to strong spatio-temporal coupling. The formation of a 
ring around the on-axis peak, a feature not observed in 
the collimated case, is shown in Fig. (b) and indicates 
steepening and breakup of the electric field. 

Note that up to now, all temporal field profiles are plot¬ 
ted on-axis, however at different positions in the radial 
dimension we expect to have more or less shock forma¬ 
tion in the temporal electric field profile. Fig. shows 
the electric field of a wavepacket, propagated 170 jum in 
diamond without any initial curvature, as a function of 
the radial coordinate. As we can clearly see, the amount 
of field steepening is directly proportional to the field 
amplitude, and therefore varies in the radial dimension. 
The center part of the beam is showing the ” top-hat” 
field shape, while as we move outwards the field steep¬ 
ening progressively weakens. At the outer part of the 
beam, the field retains it sinusoidal shape almost unal¬ 
tered. This trend should be taken into account for cases 
where the beam profile undergoes spatial reshaping, as 
in the case shown in Fig. 0(b, d). 




FIG. 9. (Golor online) (a) Broadened spectrum of the 8-/xm 
wavepacket after propagation in single-crystal diamond (cyan 
continuous line). Dashed lines: artificially truncated spectra, 
and (b) corresponding electric fields. 


V. FIELD SHOCK CHARACTERIZATION 

In this section a detailed study and characterization 
of the carrier wave shock will be conducted. As was 
stressed in our previous work [4], the main contributor 
to the formation of a carrier shock is the generation and 
co-propagation of the third harmonic wave. Here we will 
analyze and justify this claim in more detail. The carrier 
shock dispersion length scale Lcsd was introduced in 
II B| as a natural length scale of the MKP equation. We 
also demonstated that Lcsd is directly related to the 
third harmonic coherence length. In Fig. (a) we can 
see the generated higher harmonics of the 8-/xm beam 
in diamond, which is represented by the cyan continuous 
line. The corresponding field shape is shown in Fig.[^(b) 
again with the cyan continuous line. To isolate the effect 
each of the harmonics has on field shape, we artificially 
truncated the spectrum up to the 3rd harmonic, shown in 
Fig. I (a) with the red squares. The corresponding field 
shape, shown in Fig. [^(b) (red continuous line), is much 
more jagged when compared to the case where all the 
harmonics are taken into account. However, the electric 
field is still steepening, and the shock is still observable. 
In addition, the field shape is only slightly changing as 
the number of generated higher harmonics is varied from 
1 to the full spectrum. On the other hand, if we keep 
only the fundamental by truncating all higher harmon¬ 
ics (including the 3rd), the field shape remains virtually 
unmodified as shown in Fig. (b) (black dashed line). 
This clearly demonstrates that the actual steepening of 
the field is mainly coming from the third harmonic, while 
the rest of the harmonics smooth out the jagged oscilla- 
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FIG. 10. (Color online) Characteristic length scales for the 
nonlinear Kerr effect (blue dotted line), carrier shock disper¬ 
sion (black dashed line), and collapse distance (red continuous 
line) in single-crystal diamond as a function of wavelength. 
The input intensity is 1 x 10^^ W/cm^. 


tions observed when only the 3rd harmonic is taken into 
account. The oscillations on the field profile are inter¬ 
preted as temporal interference fringes between multiple 
harmonics. In the case of the full spectrum, the large 
number of harmonics decrease the period of the fringes, 
effectively smoothing out the field profile. 

The above observation validates our claim that the car¬ 
rier shock dispersion length LpsPi i t was introduced 
first in [4] and derived in detail in IIB, is in fact a very 
robust way to predict shock formation by taking into 
account the material dispersion and laser pulse parame¬ 
ters. Lcsd is a way to measure the length over which 
the carrier field shock dissipates due to dispersive effects 
of the medium. The main characteristic lengths, besides 
LcsDi that are of interest here are the nonlinear length 
scale Lnl and collapse distance Lc as defined in Ref. 


Lc = ^^ ' (24) 

ViiPinlPcrV^^ - 0.852)2 _ 0,0291 

where Pin is the total input power and Per is the critical 
power for self-focusing depending on the shape of the 
wave packet. For Gaussian spatio-temporal wave packets 
it is approximately given by Per = 3.77Ao/(87rnon2). For 
all materials in this study, both L^l and Lc are related 
to the instantaneous optical Kerr effect. 

It is now possible to directly compare the significance 
of each physical effect for a given intensity value by 
their characteristic lengths. In Fig. we can see all 
three characteristic lengths plotted for a peak intensity 
of Jo = 1 X 10^^ W/cm^ (a typical value in the fila- 
mentation regime) as a function of wavelength in single¬ 
crystal diamond. The propagation of the wave packet is 
governed by the physical effect with the shortest char¬ 
acteristic length. The general trend is that at shorter 
wavelengths carrier shock dispersion is strong and the 
generated harmonics walk-off before field steepening can 
occur. At longer wavelengths, where dispersion is much 
weaker, the walk-off of the higher harmonics is slow, and 
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FIG. 11. (Golor online) (a) Gharacteristic length surfaces 
of single-crystal diamond plotted in 3D. Black grid: Lcsd, 
color-coded surface: Lnl- (b) Minimum length scale over the 
/ — A plane. Black region: Lcsd < Lnl and blue region: 
Lnl < Lcsd indicative of shock formation. 


a carrier wave shock is expected to occur. 

Note that the carrier shock dispersion length scale is 
independent of the intensity. On the other hand, L^l 
depends on intensity, which for practical reasons makes 
a direct comparison between the two difficult due to the 
ever shifting peak intensity found in the filamentation 
regime. For this reason it is helpful to scan over a range of 
intensity values and produce a surface for each Lnl, and 
Lcsd, the latter being invariant along the intensity-axis. 
Fig. [TT]^a) shows the overlapping surfaces for Lnl (color 
coded surface) and Lcsd (black grid) for single-crystal 
diamond. Since the shortest length scale dominates, car¬ 
rier shock is expected to manifest itself in areas where the 
surface corresponding to Lnl is under the surface corre¬ 
sponding to Lcsd- To simplify the 3D surfaces shown 
in Fig. prj^a) we project them onto the I — X plane. In 
Fig.[TT[bJ we can see that the I—X plane is now separated 
in two regions, a region where carrier shock is expected 
to occur (blue) and a region where it is not (black), de¬ 
pending on which of the relevant length scales is shorter. 
The boundary between the two regions is given by the 
relation 


iC) 


cuj\k {3uj) — 3k (cj)! 
n2 


(25) 


when Lcsd = Lnl- FigjTT] provides a practical estimate 
for whether or not carrier shock formation is expected to 
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FIG. 12. (Color online) Minimum length scale over the I — X 
plane for air, xenon, argon, krypton, helium, neon and fused 
silica. Black regions: Lcsd < Lnl and blue regions: Lnl < 
Lcsd indicative of shock formation. 


happen in diamond. 

The comparison between Lcsd and L^l shown in 
Fig. 11 (b) can be done for any material with known dis¬ 
persive and nonlinear coefficients. Fig, shows the I — X 
planes, separated into the ’’shock region” {Lj^l < Lcsd) 
in blue and the ’’collapse region” {Lcsd < L^l) in black, 
for all noble gases, air, and fused silica glass. As we 
can see, each material has a distinct region where carrier 
wave shock is expected to be observed. Furthermore, 
we emphasize that these regions of carrier wave shock 
are entirely determined by the dispersive and nonlinear 
properties of the medium. Note that not all materials 
are suited for field shock generation, both neon and fused 
silica have Lcsd < L^l in the I — X region of interest 
here. Neon and fused silica represent two extremes, the 
former being very weakly nonlinear, the later being very 
dispersive. It should be stressed that Lcsd is a length 
scale corresponding to the MKP equation in a normal dis- 
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FIG. 13. Maximum of derivatives of the electric field for air 
(top) and single-crystal diamond (bottom). Propagation is 
simulated independent of a radial coordinate such that self- 
focusing and diffraction do not impact the results. 


persion regime and that it primarily quantifies how fiat 
the dispersive landscape is at a given carrier wavelength. 
For media with anamalous dispesion, such as fused silica, 
there are wavelengths where Lcsd is no longer a useful 
predictor of shock [37]. 

Since carrier wave shock formation is mainly driven by 
competition between nonlinearity and dispersion, the lat¬ 
ter being a fixed material property, it is expected that the 
amount of steepening in the electric field will vary along 
the propagation length, following the peak intensity to 
some extent. In Fig. we can see the maximum of the 
temporal derivatives of the electric field as a function of 
propagation distance for air and single-crystal diamond. 

The slight field steepness saturation after propagating 
50 cm in air, and 120 fim in diamond, occurs because the 
harmonics constantly walk-off. The field undergoes self- 
steeping a multitude of times as is indicated by the nu¬ 
merous oscillations for both materials. This happens be¬ 
cause of the dynamic interaction between harmonic gen¬ 
eration and walk-off. The oscillation period varies with 
intensity and material dispersion but is always observed. 
Similar results have been seen in all materials studied 
here (data not shown) as well an in the preceding work 
|4], indicating that the phenomenon quite robust. 


VI. ROBUSTNESS AND UNIVERSALITY 

A very important observation is that the carrier shock 
formation is unaffected by the pulse duration since nei¬ 
ther Lcsd nor Lnl depend on it. In Fig. we can 
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FIG. 14. (Color online) Electric field of a Tp = 200 fs, Ao = 8 
fim pulse, propagated in single-crystal diamond for 170 /xm. 


see the electric field of a 200 fs pulse undergoing self- 
steepening in diamond after 170 /xm of propagation. The 
carrier wave shock forms over multiple oscillations of the 
electric field. Interestingly, the reshaping of the electric 
field is not constant over the whole pulse, but rather the 
steepening follows the field amplitude. Oscillations at 
the front and back of the pulse show very little steepen¬ 
ing, while as one moves closer to the central time slice, 
’’shark-fin” like features become evident. At the center 
of the pulse the electric field is strongly reshaped into a 
’’top-hat”, effectively having gone through all stages of 
the carrier shock formation. As stated before, the carrier 
shock formation is intensity (or else field amplitude) re¬ 
lated and is acting on each point of the field separately 
through the co-propagation of the generated harmonics. 
Therefore, the overall pulse duration is of little impor¬ 
tance concerning carrier shock formation. However, the 
use of longer duration high power pulses will be limited 
by optical damage in solid state media like single-crystal 
diamond. 

As mentioned earlier, we will investigate the role of vi¬ 
brational Raman in the case of single-crystal diamond us¬ 
ing the model described in section III. In order to isolate 
the effect of vibrational Raman optical field ionization, 
plasma is switched off. Simulations were conducted using 
the UPPE model. Conclusions drawn here should apply 
to all materials with a clearly defined Raman response. In 
Fig. we can see the effect of vibrational Raman on the 
carrier shock formation of an 8 /xm - 40 fs pulse in single¬ 
crystal diamond. Although the frequency and time scale 
are known [34], the actual relative strength (in compari¬ 
son to the instantaneous Kerr effect) at the given wave¬ 
length is not. Therefore the relative strength is varied 
within a reasonable range from 0% (Fig. [T^a)), to 25% 
(Fig. [^b)), up to 50% (Fig. [^c)). In Figs, [l^a)-(c) 
we observe that as the Raman strength is increased the 
amount of carrier shock is decreasing. However, the effect 




FIG. 15. (Golor online) Electric field profiles for various Ra¬ 
man strengths after propagation of z = 80 /xm in single-crystal 
diamond. Raman contribution of (a) 0% , (b) 25%, (c) 50%, 
and (d) 50% with double the intensity value used in (a)-(c). 


of Raman on the field shock formation is mostly neutral 
rather than destructive. This becomes evident when we 
increase the intensity by a factor of two in the case of 
50% Raman strength, which effectively reverses the loss 
of field shock due to the inclusion of Raman. There¬ 
fore, simply increasing the intensity of the pulse should 
in principle still lead to significant shock formation, given 
that intensity is kept under the damage threshold of the 
medium. This would mean that given precise knowledge 
of the Raman relative strength, a correction has to be 
applied to the shock formation region in the I-A plane in 
Fig. and Fig. This would result in a shift of the 
collapse region (black) upwards, effectively shrinking the 
shock-formation region (blue) by an amount depending 
on the strength of the Raman response. 

Lastly, we will study the effect of high nonlinear losses 
and plasma, including plasma defocusing and absorption, 
on carrier shock formation. Even though ionization is 
properly modeled as shown in section HB. and Fig.[^ an 
accurate and well established ionization model for single¬ 
crystal diamond at mid-IR wavelengths is not available in 
the literature. Therefore, we investigated and compared 
carrier wave shock formation for two extreme cases: us¬ 
ing artificially high ionization coefficients and switching 
ionization off altogether. 

In Fig. we observe the effect of artificially increas¬ 
ing ionization on the carrier shock formation of our 8/xm 
- 40 fs wavepacket in single-crystal diamond. As we can 
see, the ’’shark-fin” and ’’top-hat” shaped electric fields 
are significantly reduced in amplitude (Fig. 16 (c), (d)) 
when compared to the case without ionization shown in 
Fig. 16 (a), (b). However, even with artificially high ion¬ 
ization coefficients, field steepening is still observed. Ion- 
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FIG. 16. (Color online) Electric field profiles of an 8 /xm - 40 fs 
pulse in single-crystal diamond, without (first row) and with 
ionization (second row). Black line: initial field shape, red 
lines: reshaped electric field after z = 80 /xm (first column) 
and = 200 /xm of propagation (second column). 


ization losses tend to reduce field amplitudes, and conse¬ 
quently they reduce the amount of field steepening. On 
the other hand, plasma generation has been associated 
with a strong blue-shift in the spectrum and supercontin¬ 
uum generation [2] , which could in principle contribute to 
field shock formation in a similar way as higher harmonic 
generation does. The effect of plasma on field shock for¬ 
mation is proving to be a complex process, depending on 
setup parameters and propagation dynamics. We have to 
stress again that the use of shorter wavelengths will lead 
to an increase of ionization, which along with stronger 
dispersion makes carrier shock formation extremely diffi¬ 
cult to observe. On the opposite side, moving towards 
longer wavelengths favors field shock formation, since 
ionization is weaker and Lcsd becomes longer. 
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FIG. 17. (Golor online) Gomparison of carrier shock forma¬ 
tion between GMKP (black dashed lines) and UPPE (red 
continuous lines), for an 8 /xm - 40 fs wavepacket in single¬ 
crystal diamond, (a) Generated ’’shark-fin” electric field for 
UPPE and GMKP. (b) Gorresponding spectral intensities, (c) 
Group velocity dispersion as modeled by GMKP and UPPE 
for single-crystal diamond vs wavelength A, and (d) vs nor¬ 
malized frequency. 


model is fully justified. In Figs. (c) and (d), the dis¬ 
persion relation of diamond is shown for both models, 
visualized in wavelength and normalized frequency re¬ 
spectively. As we can see there is a slight discrepancy 
between the GMKP dispersion and the actual Sellmeier 
dispersion relation used in the UPPE model, especially 
at high frequencies. This should be taken into account 
when conducting numerical experiments at shorter wave¬ 
lengths. For shorter wavelengths, the GMKP model can 
poorly approximate material dispersion and hence sig¬ 
nificant differences between the UPPE and GMKP mod¬ 
els may result. This means that the pulse’s wavelength 
can determine the numerical model of choice, UPPE for 
shorter wavelengths where dispersion is accurately mod¬ 
eled and GMKP at longer wavelengths since it is the 
canonical equation describing pulse propagation in that 
regime [4]. 


VII. COMPARISON BETWEEN GMKP AND 
UPPE 

To validate our result we compared the results of 
the GMKP model with the full UPPE field propagator 
throughout this work. In Fig. we can see minor dis¬ 
crepancies between the two models in the case of our 8 /xm 
- 40 fs wavepacket propagating in single-crystal diamond. 
For simplicity, no vibrational Raman or optical field ion¬ 
ization were taken into account. In Fig.[^(a) we can see 
the predicted carrier wave shock using both the GMKP 
(blue lines) and UPPE (black dashed lines) models. The 
corresponding spectral intensities are shown in Fig. 

(b). The field profiles are essentially indistinguishable, 
showing that at long wavelengths the use of the GMKP 


VIII. CONCLUSION 

We have presented a broad and detailed numerical 
study of carrier wave shock formation in the weakly dis¬ 
persive - highly nonlinear regime in air, noble gases and 
single-crystal diamond. ’’Shark-fin” and ”top-hat” elec¬ 
tric fields are predicted to form spontaneously at long 
wavelengths for a variety of input powers and pulse du¬ 
rations. In addition, the field shock formation is proven 
to be very robust since it is occurring despite ionization 
losses, plasma generation, vibrational Raman and two- 
phonon absorption. Carrier shock is predicted to be un¬ 
avoidable in most high intensity, long wavelength pulse 
propagation settings. The results of the GMKP model 
were compared with the ones from the full field UPPE 
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model, and were found to be in good agreement with 
each other. Results presented here are expected to have 
significant impact on multiple disciplines in the field of 
nonlinear optics, especially attosecond pulse generation 
na, waveform synthesis m, and pulse compression [12] . 
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